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Abstract 

We present the first implementation in the 7) plane of the generator coordinate method with 
full triaxial angular momentum and particle number projected wave functions using the Gogny 
force. Technical details about the performance of the method and the convergence of the results 
both in the symmetry restoration and the configuration mixing parts are discussed in detail. We 
apply the method to the study of ^^Mg, the calculated energies of excited states as well as the 
transition probabilities are compared to the available experimental data showing a good overall 
agreement. In addition, we present the RVAMPIR approach which provides a good description of 
the ground and gamma bands in the absence of strong mixing. 

PACS numbers: 21.60.Jz, 21.10.Re, 21.60.Ev, 27.60.+j 
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I. INTRODUCTION 



Self-consistent mean field methods with effective phenomenological interactions and their 
extensions beyond mean field provide the appropriate theoretical tools for describing many 
phenomena along the whole chart of nuclides, from light to medium, heavy and superheavy 
nuclei in or far away from the stability valley [1]. On the one hand, the success of these 
methods is related to the high quality of the phenomenological effective interactions used 
-Skyrme, Gogny or Relativistic Mean Field (RMF). On the other hand, the mean field 
method allows the inclusion of many correlations within a very simple intrinsic product 
wave function. Hence, bulk properties such as masses and radii are very well described at 
the mean field level. However, in some cases this picture fails to take into account important 
correlations and methods beyond the mean field approach have to be applied. Furthermore, 
because the mean field is defined in the intrinsic frame it is mandatory to go beyond this 
approximation to evaluate excitation energies or transition probabilities in the laboratory 
system. 

There are several methods to incorporate the correlations missing at the mean field level. 
Normally, the intrinsic wave functions are allowed to break relevant symmetries of the sys- 
tem, for example, particle number, rotational and translational invariance, parity, time- 
reversal, etc. to enlarge the variational space and incorporate, for instance, deformation or 
superfluidity in the mean field picture. This leads to a degeneracy of the wave functions 
rotated in the gauge space associated to the broken symmetry. An appropriated superposi- 
tion of these wave functions provides a symmetry conserving many-body wave function and 
an additional lowering of the energy of the system. In this way, using projection techniques 
[2], many correlations are obtained by restoring some or all of these symmetries . Further- 
more, the mixing of different mean field configurations within the general framework of the 
Generator Coordinate Method (GCM) [2] allows the inclusion of quantum fluctuations along 
some relevant collective variables such as the multipole moments. 

Most of the currently used beyond mean field calculations with effective forces include two 
symmetry restoration, i.e., particle number (PN) and angular momentum projection (AMP) 
and configuration mixing along the axial quadrupole deformation [1, 3, 4]. This approach 
(axial GCM-PNAMP) has been successfully applied to study many phenomena like, for ex- 
ample, the appearance or degradation of shell closures in neutron rich nuclei [3, 5-7], shape 



2 



coexistence in proton rich Kr [8] or Pb [9, 10] isotopes or shape transitions in the A ~ 150 
region [11, 12]. However, the intrinsic wave functions used there were restricted to have 
axial symmetry, with K = 0, because this assumption simphfies considerably the angular 
momentum projection and lightens the computational burden significantly. This restriction 
is one of the major drawbacks of the method because it limits its applicability to systems 
where triaxiality does not play an important role. However, many exciting experimental and 
theoretical phenomena are closely related to the triaxial degree of freedom, for instance: the 
presence of 7-bands in the low lying energy spectra and 7-softness, shape coexistence and 
shape transitions in transitional regions [13, 15-20]; the lowering of fission barriers along the 
triaxial path [21-23]; the influence of triaxial deformation in the ground state for the mass 
models [24, 25]; triaxiality at high spin [26-28]; the observation of ii"-bands and isomeric 
states in the Os region [29-31]; or some other exotic excitation modes such as wobbling 
motion and chiral bands [32-34]. 

From the theoretical point of view some approaches beyond mean field have been proposed 
to study the triaxial effects. In particular, one of the most widely used is the collective 
Hamiltonian [2] given in different versions depending on the underlying nucleon-nucleon 
interaction used to define the collective potential, namely Pairing-plus-Quadrupole [35], In- 
teracting Boson Model [36], Nilsson Woods-Saxon [26], Gogny [37-39] or RMF [40]. This 
model has been applied successfully to describe some of the experimental features men- 
tioned above. However, the collective Hamiltonian can be understood as a gaussian overlap 
approach (GOA) of the triaxial GCM and this description should be improved including 
properly the effects of the symmetry restoration and the full configuration mixing without 
any GOA approximation. 

In the past, exact angular momentum projection with triaxial intrinsic wave functions 
without GCM have been carried out only for schematic forces and/or reduced configura- 
tion spaces. Examples are the projection of BCS [41] or Cranked Hartree-Fock-Bogoliubov 
(CHFB) states [42] with the Pairing-plus-Quadrupole interaction, the projection of Cranked 
Hartree-Fock (CHF) states (no pairing) with schematic [43] and full Skyrme interactions [44] 
or angular momentum projection before variation with particle number and parity restora- 
tion in limited shell model spaces [45, 46] . 

However, recent improvement of the computational capabilities enabled the first implemen- 
tations of the angular momentum projection of triaxial intrinsic wave functions in the whole 
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(/5,7) plane with effective forces. In particular, Bender and Heenen reported GCM cal- 
culations with particle number and triaxial angular momentum projection (PNAMP) with 
the Skyrme SLy4 interaction [47]. In this work, the intrinsic wave functions were found by 
solving the Lipkin-Nogami (LN) equations. On the other hand, Yao et al presented the 
implementation of the triaxial angular momentum projection [48] and the extension to the 
GCM [49] for the Relativistic Mean Field (RMF) framework. In the latter work no particle 
number projection has been performed and the mean field states are found by solving the 
RMF+BCS instead of the full HFB or LN equations. These approximations could lead to 
a poor description of important pairing correlations, especially in the weak pairing regime 
where even spurious phase transitions appear [3, 50]. 

In this paper we present the first implementation of the Generator Coordinate Method with 
Particle Number and Angular Momentum Projected (GCM-PNAMP) triaxial HFB wave 
functions with the finite range density dependent Gogny force [51]. The finite range of the 
Gogny force provides excellent pairing properties and is often used as a benchmark in this 
respect. Furthermore it is able to provide at the same time both good global as well as spec- 
troscopic properties [52, 53]. The intrinsic HFB states are found by solving the Variation 
After Particle Number Projection (VAP-PN) equations [54]. This fact constitutes the main 
methodological difference with respect to the calculations reported in Ref. [47]. This is a 
very important difference because VAP-PN allows the inclusion of the pairing correlations 
in a very efficient way yielding a significant improvement of the final results with respect to 
other approaches [3, 54]. 

In nuclei without strong mixing the so called Variation After Mean-field Projection In Real- 
istic model spaces (VAMPIR) [45] approach has been very successful. In this approach only 
one HFB wave function is considered which is determined by minimization of the projected 
energy, i.e. the VAP approach for both the AM and the PN projections. A full VAMPIR 
approach with the Gogny force and large configuration spaces is not feasible yet. Instead we 
use an approximation to it, which we call RVAMPIR and which as we shall see, in the case 
nucleus studied in this article, provides very reasonable results for the ground and 7 bands 
with much less effort than in the GCM approach. 

The paper is organized as follows. In Sec. II we will give an overview of the theoretical 
framework. Then, we will focus our analysis on the nucleus ^^Mg which has been studied 
as a test case in earlier implementations of the GCM-PNAMP method with Skyrme and 
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Relativistic interactions. In particular, in Sec. Ill we will show a standard axially symmetric 
calculation which allows to make an educated guess for some relevant parameters needed in 
the full calculation such as the number of major oscillator shells or the relevant deforma- 
tions ranges. In Sec. IV we will analyze in detail the simpler PNAMP method, studying the 
convergence of the integrals in the Euler angles, giving some consistency requirements and 
showing the role of having an adequate mesh in the (/3,7) plane. In Sec. V we will show 
the final results for the calculated spectrum and B(E2) transitions strengths of ^^Mg and a 
comparison with experimental data. Finally, a brief summary and outlook on future work 
will be addressed in Sec. VI. 

II. THEORETICAL FRAMEWORK 

A. Generator Coordinate Method with Particle Number and Angular Momentum 
Projected states (GCM-PNAMP) 

In the present approach, the final many-body wave functions that describe the different 
states of an even-even nucleus with Z{N) number of protons (neutrons) are written as: 

|/M; NZa) = f^^^^'''\IMK; NZ; (1) 

KM 

where (/9,7) are quadrupole deformation parameters (see below), a = 1,2,... labels the 
levels for a given value of the angular momentum / and M, K are the projections of / 
on the laboratory and intrinsic z— axes respectively. The coefficients fll^^'" of the linear 
combination are found by minimizing the energy within the non-orthogonal set of wave 
functions {\I M K ; N Z . These states are obtained by projecting the intrinsic mean- 
field states |$(/3,7)) onto good particle number and angular momentum: 

\iMK-Nz-p^) = J I?[*^(fi)J^(^])p^p^|<|.(/3,7))rf^] (2) 

with = 2^ Jq^ e^'f'^^~^^dip the neutron number projector ( the associated gauge angle 
and protons the proton number projector), i?(f2) and P{j^(i7) are the rotation operator 
and the Wigner matrices [55] in the Euler angles VL = (a, b, c) [69], respectively. In principle, 
the ranges for these angles are (0 < a < 27t, < b < tt,0 < c < 2tt). However, for intrinsic 
Hartree-Fock-Bogoliubov (HFB) states (|$(/3,7))) which are symmetric under time-reversal 
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and simplex symmetries, the intervals for both gauge and Euler angles can be reduced to 

(0 < < 7r/2) and (0 < a < 7r/2, < 6 < 7r/2, < c < vr), respectively [47]. 

The wave functions (Eq. 2) are eigenstates of the particle number and angular momentum 

operators: 

N\IMK;NZ;l3-f) = N\IMK; NZ; jS-f), (3) 
Z\IMK;NZ;(3-f) = Z\IMK; NZ; (4) 
P\IMK;NZ;(3-f) = h^I{I + 1)\IMK; NZ; l3-f) (5) 
i^\IMK;NZ;l3-f) = hM\IMK; NZ; jS-f), (6) 
h\IMK]NZ][5-f) = hK\IMK;NZ;l3-f) (7) 



The intrinsic HFB states (|$(/3,7))) are obtained by minimizing the particle-number pro- 
jected energy functional -E^'^ $(/3,7) (variation after projection, VAP) [54]. This is one of 
the most relevant parts in the calculation because the quality of the result largely depends 
on the structure of the intrinsic HFB-type wave functions used. In contrast to other methods 
like plain HFB or Projected Lipkin-Nogami (PLN), the VAP-PN performs the restoration 
of the particle number symmetry in an optimal way, including pairing correlations both in 
the weak and strong pairing regimes [3]. This is especially relevant in GCM-like theories 
where a large grid of (/3,7) points is needed. The strength of the pairing correlations has 
a strong dependence on the single particle level density and the latter one itself with the 
deformation parameters. This implies that a strongly (/3,7) dependent oscillating pairing 
regime appears in the calculations and consequently theories like plain HFB (BCS) or PLN 
(LN) are unable to cope with this challenge providing wave functions of oscillating goodness. 
Only a VAP-PN approach warrants high quality solutions independently of the (/3, 7) values. 
Dealing with effective forces like Skyrme, Relativistic and Gogny, a natural separation of the 
interaction into the two-body Hamiltonian on the one hand and the density- dependent 
part, £:^^[$] on the other emerges. In our case, we are using the Gogny DIS interaction 
[51] and H2h corresponds to the kinetic energy (the two-body part from the center of mass 
correction included) plus the spin-orbit. Coulomb and the finite range central potentials. 
In the calculations, all direct, exchange and pairing terms are included [56]. The VAP-PN 
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principle, provides 

[<l(/3,7) 

where: 



_ = (8) 



^""'^[•^l = JLr>zl + ^^^^t*^] - V(<^l4o|$) - A,,,($|42|$) (9) 
($|P^P^|$) 

In a beyond mean field method, and in particular for the particle number projection, we 
need a reasonable prescription for the spatial density, which we shall call Pint(^)! that enters 
in the density dependent term of the interaction. In this work, assuming the 

phenomenological nature of these interactions and considering that the restoration of the 
particle number symmetry is performed not in the coordinate but in the gauge space, we 
have chosen the number projected spatial density prescription that has proven to be free 
of divergences [54] and to give very good results for describing many phenomena along the 
nuclear chart: 



with p(r) = / df'6{f — r'). As shown in [54] for the PNP and in [57] for the Lipkin-Nogami 
approach, the use of the projected density or the so-called mixed prescription (in the case 
when the latter is free of potential divergences) provide very similar results. Furthermore, 
we see in Eq. 9 that the minimization is performed under constraints on the quadrupole de- 
formation operators Q2^- The Lagrange multipliers Xq^^ ensure that the following conditions 
are fulfilled in the intrinsic state: 

Ko ('^'IQ2o|'^) = q2o 

A,,, ^($|Q22|*)=g22 (11) 

In addition, the deformation parameters (/3,7) are directly related to (520, ^722) by: 

_/3cos7 _/3sin7 fV An 

being ro = 1.2 fm and A the mass number. These constraints allow to explore the (/3,7) 
plane to generate the wave functions to be used in the configuration mixing calculations. 
We now describe the Generator Coordinate Method (GCM) to obtain the final spectrum 

■I;NZ 
if/37 

if/37 

E(njI]NZ rpI;NZ;o- \rI\NZ \ rI;NZ;a _ r. /-, o\ 

iC'/3'7' 



^^i;NZ;a-^ and the coefficients fK^^^'" given in Eq. 1. Minimization of the energy with respect 
to the coefficients /^aL^''^ leads to the Hill- Wheeler- Griffin (HWG) equation 



which has to be solved for each value of the angular momentum. The GCM norm- and 
energy-overlaps have been defined as: 

A/XSw = {IMK-NZ;f3j\IMK'-NZ-f3W) 

^ W'/3'y = (^^^; HH2i.\IMK'; NZ; /3V) + ^Ld"" '"^^ [HP, 7), Y)] (14) 

In the last expression, we have separated again the energy overlap in the contribution of 
the pure Hamiltonian part of the interaction and the density-dependent term. In the lat- 
ter, we have used the particle number projected spatial density combined with the mixed 
prescription for the angular momentum projection and GCM part, namely: 

,-(ar-).«<5i^^. (15) 

This prescription is suitable for dealing with the restoration of broken symmetries in the 
coordinate space such as the rotational invariance or the spatial parity. 
Once we have calculated the corresponding GCM overlaps, the next step consists in solving 
the HWG equations (Eq. 13). To cope with the problem of the linear dependence one first 
introduces a orthonormal basis defined by the eigenvalues n^' and eigenvectors u^^^.^ of 
the norm overlap: 



K'I3'Y 

This orthonormal basis is known as the natural basis and for Ux^^ values such that 
n^'^ jn\^^ > (, the natural states are defined by: 

I;NZ 

|^7M;A.z^ = E ^7^|/Mir;iVZ;/37). (17) 

Obviously, a cutoff ( in the value of the norm eigenvalues has to be introduced in order 
to avoid linear dependences [48]. Then, the HWG equation is transformed into a normal 
eigenvalue problem: 



A' 

"^A 



From the coefficients Ci'^^''^ we can define the so-called collective wave functions 



7) that account for the probability density, normalized to 1, of finding the state 
(/, cr) with given deformation parameters (/3,7): 

7) = E G^'^^v^^^^A = E Fl^''"''iP, 7). (19) 

A,K K 



we have also introduced 7) that account for the probabihty density of finding the 

state (/, 0") with given values of K and deformation parameters (/3,7). 
Furthermore, the expectation value of a generic operator O is given by 

I;NZ* I;NZ 

^I;NZ;.^Y: E G'^''^--^*^^^{w\d\w')''''^^ (20) 
A;A' Ki3r,K'(5'Y \J n'^^^ y Upi!^^ 

with {w\d\w') = {IMK;NZ;(3'j\d\IMK';NZ;(3'y). This expression can be generalized 
to account for transitions associated to the tensorial operator T12: 

h-NZ* I2;NZ 

tiha, ^ I,a,) = J: E Gr^-'^^^^^^iu^m^llu.'^^^ (21) 
A;A' KPr,K'i3'Y ^Jn^' \J nl; 

where (ti7i| IT121 1^2) — (hK; NZ; P'y\\Ti2\\l2K'; NZ; 13'^') stands for the reduced matrix ele- 
ment calculated according to the Wigner-Eckart theorem [2, 55]. Detailed expressions for cal- 
culating these reduced matrix elements for B(E2) transitions and spectroscopic quadrupole 
moments within this framework can be found elsewhere [5, 47, 49]. 

B. Simpler approaches: Particle Number and Angular Momentum Projection 
(PNAMP) and the RVAMPIR approximation 

The expressions given above constitute the most general framework that we are using for 
solving the nuclear many body problem. Nevertheless, there are some limiting cases with a 
relevant physical meaning that can be deduced in a straightforward manner from them. The 
first one is the particle number projection (PNP) that has been discussed above (Eq. 9). 
The second approach is the particle number and angular momentum projection (PNAMP) 
of a single point in the (/3,7) plane. Here, the wave function is of the form of Eq. 1 but 
without the mixing in the deformation parameters: 

\IM-NZ-u-p,^) = Y.h'f^^''{M)\IMK-NZ-p^), (22) 

K 

where the label v stands for the (2/ + 1) different states that can be obtained with the 
angular momentum projection. However, due to the time reversal and simplex symmetries 
imposed on the intrinsic wave functions, this number is reduced to (//2 + 1) and ((/ — l)/2) 
states for even and odd values of J, respectively. Moreover, if we furthermore have axial 
symmetry, only one state can be obtained and only for even values of /. 
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The coefficients /i}^ '""(/^t) and the PNAMP energies E^'^^''^{P,'j) are found by solving the 
simphfied version of the HWG equation (see Eq. 13): 

E - E^-^'^'-'^P. iW'k^,%,,) hT'-{P. 7) = 0. (23) 

K 

The remaining expressions used to solve the HWG equations are simplified in the same 
manner, i.e., removing the sum over (/3,7) from the equations and evaluating only the 
diagonal part. In addition, the collective wave functions F^'^^''^(/3, 7) (Eq. 19), which in 
analogy we shall call ifj^^^'^(/3, 7), now give the spectral distribution in the K space of the 
corresponding PNAMP state. 

A full variation of the HFB wave function in the VAP approach, in the spirit of VAMPIR, 
for the PN and the AM with large configuration spaces and the Gogny interaction is not yet 
feasible. However we can use an approximation to VAMPIR, which we shall call from now 
on RVAMPIR, in which the PN is handled in the VAP approach and the AM in a Restricted 
VAP (RVAP) one. The RVAP approximation has been thoroughly studied in [58]-[59]. 
In the VAP method the whole Hilbert space associated with the HFB transformation is 
scanned in the variational procedure. In the RVAP method, however, only a restricted 
variational space of highly correlated wave-functions is allowed in the minimization process. 
Monopole (pairing) and quadrupole (/3 and 7) correlations are believed to be the most 
relevant degrees of freedom of atomic nuclei and are related to the particle number and the 
angular momentum symmetries, respectively. Since we are considering the PN symmetry in 
the VAP theory it seems reasonable in our case to consider the restricted Hilbert space to 
contain a whole set of quadrupole deformed wave-functions |$(/3,7)) which parametrically 
depend on (/3,7). This procedure is justified by theoretical arguments [2] which establish 
that a VAP approach is needed for systems with weakly broken symmetries, like in the 
PN case where only a few Cooper pairs participate, but it can be approximated in case 
of strongly broken symmetries, such as deformation, where a large number of nucleons 
participate. Concerning the differences of this approximation as compared to VAMPIR it is 
clear [58] that if, besides of considering the quadrupole moments Q20 and Q22 in Eq- 9, we 
will include higher multipole moments Qlm to increase the variational space, our solution 
would get very close to the one of the genuine VAMPIR. With respect of the quality of our 
approach (again with respect to the full VAMPIR) we expect that in general it will be very 
similar and only in very soft nuclei, where higher modes (hexadecupole for example) are 
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very relevant, differences may arise. But for very soft nuclei we have to question also the full 
VAMPIR since a GCM-like approach will be more appropriate. That means, RVAMPIR is 
not as "restricted" as its name might imply. 

Specifically the basic RVAMPIR approach consist of the following steps: 

A. - At each {(3^,%) value of a given set of points in the (/3,7) plane the following items 

are performed: 

Al.- Solve the VAP-PN equations, Eqs.8-9, to determine the /3 — 7 constrained HFB 
wave function |$(/3,7)). 

A2.- Carry out simultaneous particle number and angular momentum projection on 
the wave function |$(/3,7)), what we have called \IMK; NZ; P'j) , see Eq. 2, to 
form the linear combination of the state \IM; NZ; u; fS'y) in Eq. 22. 

A3.- Solve the HWG equation, Eq.23, for different angular momenta. 

B. - For each value of the angular momentum sort out the energies E^'^^'^IP, 7) of Eq. 23 

and find out the point (/34in,7min) providing the energy minimum £^min^'''(/34in' T^in) 

C- The solutions of the HWG equation at the points {(3^^^, 7min) provide //2 + l((/ — 1)/2) 
states for even (odd) / values, which allow to build a partial spectrum and to calculate 
the transition probabilities among the different states or any other observable. 

One has to notice that all RVAMPIR states are orthogonal, those with different AM 
in an obvious way and those with the same AM because they are solution of the same 
eigenvalue equation. 

In the following sections we will give some examples of the convergence, consistency and 
performance of the methods described above. All the many body intrinsic wave functions 
and operators have been expanded in a cartesian harmonic oscillator single particle basis 
closed under rotations [60]. In particular, the rotation operator R{Q) has been evaluated 
following the expressions given in Ref. [61] and the Neergard method [62] has been used 
in the calculation of the norm overlaps in order to determine the correct sign of the 
Onishi formula [63-65]. The overlaps of a generic operator have been calculated using the 
generalized Wick theorem [64]. 
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III. AXIAL CALCULATIONS FOR ^^MG 



Due to the huge computational cost of the full triaxial calculation, it is important to 
study first the axial case (with K = 0) in order to fix some relevant quantities. The most 
important ones are the region of /3 deformation to be included in the calculation and the 
number of major oscillator shells in which the mean field wave functions are expanded. The 
computational effort depends critically on these quantities and it is important to ensure the 
convergence of the results, at least in the axial case, to have a reasonable choice which then 
later allows to perform the full triaxial calculation. 

The main advantage of considering only axial symmetric {K = 0) intrinsic wave functions 
|$(/3,7 = 0°, 180°)) = |$(/3)) is that the integration over the Euler angles (a,c) can be 
done analytically and this fact reduces drastically the computational time. The simplified 
expressions of the axial GCM-PNAMP method can be found in detail in Ref. [5]. We first 
analyze the results obtained for the nucleus ^^Mg using Ng^^s = 7 oscillator shells and 
Npoints = 31 intrinsic wave functions distributed in the interval (—1.5 < /3 < 1.5) with posi- 
tive and negative values of /3 corresponding to prolate 7 = 0° and oblate 7 = 180° shapes. 
The integration over the gauge angle if for the particle number projection part has been 
performed using the Fomenko expansion [66] while for the integration over the Euler angle b 
a Gaussian-Legendre quadrature has been used. We have chosen Npom = 9 and Nf, = 16 as 
the number of integration points for the particle number and the angular momentum parts 
of the projection, respectively. With these assumptions the expectation values for the N, Z, 
iV^, Z"^ and P operators differ by less than 10"*^ from the corresponding eigenvalues. In Fig. 
1(a) we plot the potential energy surfaces (PES) along the (3 direction for the VAP-PN and 
PNAMP approaches. The VAP-PN curve shows two differentiated minima separated by a 
barrier of ~ 7.7 MeV, the lowest one at prolate deformation {/3 = 0.5) and the other one in 
the oblate part (/3 = —0.2) . These minima are shifted towards larger deformations when the 
angular momentum projection is performed. In particular, a well defined prolate minimum 
appears at /3 = 0.6 for J = 0, 2, 4, 6, 8 showing that a rotational band will develop from this 
intrinsic state. For the ground state the gain in correlation energy due to the restoration of 
the rotational symmetry amounts to ~ 3.9 MeV. In the oblate part, we observe a minimum 
at P = —0.4 for / = 0,2 while for higher values of the angular momentum the minimum 
vanishes. Here, the energy difference between the VAP-PN and the / = oblate minima is 



12 



0.4 




-1.5 -1 -0.5 0.5 1 1.5 
P 



FIG. 1: (a) Potential Energy Surfaces (PES) along the (3 deformation for particle number projection 
and particle number and angular momentum projection (^^Mg). The bullets correspond to the 
excitation energies for the different GCM levels (/, a) with their positions at = /3|F^'^(/3)|2. 
The energy is normalized to the GCM ground state energy {Qi). (b) Collective wave functions 
for the (cr = 1) GCM levels. The values of the ordinate axis is displaced by 0.05 with increasing 
angular momentum, (c) Same as (b) but for the (cr = 2) GCM levels and 2^ state. Positive and 
negative values of /? correspond to prolate (7 = 0°) and oblate (7 = 180°) shapes, respectively. 

~ 4.6 MeV. 

The next step in the calculation is the configuration mixing of the PNAMP states. Hence, 
once the HWG equations are solved, we select as the final solutions those that belong to a 
plateau in the energy as a function of the number of states in the natural basis (Eq. 17) and 
fulfill the ort honor mality condition. To avoid duplications a detailed discussion on these 
issues is postponed to the triaxial case. The resulting GCM-PNAMP energies are also rep- 
resented in Fig. 1(a), while the corresponding collective wave functions (Eq. 19) are plotted 



13 



18 - 

14 r 

_12 - 2* 2* 

?io^ . °* 

S ^ 6* 

^ 8 - 
LU 

6 - 



4 1 - 4' 

2 r 2+ 2* 

0^ °* °* 

Axial 7 shells Axial 1 1 shells 



FIG. 2: Axial PNAMP-GCM excitation spectra of ^^Mg obtained considering 7 shells (left) and 
11 shells (right) in the calculations. 

in Fig. 1(b) for a = 1 and in Fig. 1(c) for a = 2 and 2^. In these figures we can see 
that the a = 1 states are members of a rotational band, with most of the intensity of the 
collective wave functions concentrated around /3 = 0.6. This deformation corresponds to 
the location of the prolate minima of the different potential wells. The ground state Oi' also 
has a small mixing with the oblate minimum at /3 = —0.5. The situation is rather different 
for the a = 2 states. The second O"*" state is a mixing of oblate and prolate configurations, 
while wave function of the 2^ state peaks in the oblate minimum of the corresponding PES 
and the state could be considered as a vibration built on the 7 = 4 prolate well with 
a small contribution of slightly oblate states. In the state the prolate deformations are 
again favored. Remembering that the purpose of this axial calculation is to determine the 
range of values of /3 needed to obtain converged results in the low-lying energy spectrum, 
we observe in Fig. l(b)-(c) that all the collective wave functions studied here drop to zero 
at the boundaries. A smaller interval, however, could not be sufficient for describing cor- 
rectly the collective states. In addition, we have checked the convergence of the results as a 
function of the number of the points included in the GCM-PNAMP. Increasing this number 
to Npoints = 61 still yields very similar results for the PES, GCM-PNAMP energies and the 
collective wave functions as compared to the ones obtained for Npoints = 31. 
Finally, in order to test the convergence with the number of oscillator shells, we have per- 
formed a calculation with A^s/ie//s = 11 and Npoints = 31. It is noteworthy that for a triaxial 
calculation, the computational time for Nsheiis = 11 is ~30 times larger than the one used 
for Ngheiis = 7. Although this fact complicates the applicability of this method for heavy 
nuclei, for lighter systems the calculation with a smaller number of oscillator shells could 
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still be sufficient. This is the case for ^ Mg, where the PES and the collective wave functions 
calculated with Nsheiis = 11 (not shown) are very similar to the Ngh^s = 7 results. In Fig. 2 
we compare the spectra obtained in the two calculations and observe a relative error of less 
than 10% for all the levels. While the members of the cr = 1 bands almost match each other, 
small differences are found in the cr = 2 band. This comparison justifies that all further 
calculations are performed withiVs/jez/s = 7. 



IV. CONVERGENCE AND CONSISTENCY OF THE TRIAXIAL PNAMP 

In this section we will study some aspects of the simultaneous particle number and an- 
gular momentum projection with triaxial shapes. Firstly, it is important to note that the 
parametrization of the quadrupole deformation in terms of (/3,7) variables gives a triple de- 
generacy in the range 0° < 7 < 360° if we consider time-reversal conserving wave functions 
[2]. This degeneracy corresponds to the three possible orientations of the intrinsic axis /s 
with respect to the z— axis (see Fig. 3). Therefore, the interval 0° < 7 < 60° covers all the 
possible quadrupole deformations. However, we can take advantage of this symmetry first 
to improve the convergence of the integral in the Euler angles that must be carried out in 
the PNAMP calculation (Eq. 2) and second to perform consistency checks of the results. 
We now study the convergence of the integral in the Euler angles with respect to the number 
of integration points in = (a,6, c). We have considered the symmetries of the intrinsic 
wave function reducing the integration interval to (0 < a < 7r/2,0 < b < 7r/2,0 < c < tt) 
(see Refs. [42, 47, 48]) and we have used Gaussian-Legendre quadratures for the numerical 
integration. As in the axial case, the number of integration points for the particle number 
projection is kept to Npom = 9, which is sufficient to get eigenstates of the particle number 
operators. Naturally, the best candidate to check the convergence of the angular momen- 
tum projection is the expectation value of the total angular momentum operator P that, 
considering Eq. 5, must be: 

Iv'j^Ki^){<i>\R{n)p^p^\<i>)dn ^ ^ 

The convergence in the number of integration points depends on three factors, namely the 
orientation of the intrinsic axes, the values of (J, K) and the deformation /3. Let us start 
with the two latter factors. In Fig. 4 we plot the mean value of the total angular momentum 
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FIG. 3: (Color online) Orientations of the intrinsic deformation as a function of the 7 parameter. 
7 = 0°, 120°, 240° and 7 = 60°, 180°, 300° correspond to axial symmetric prolate and oblate shapes, 
respectively. 

operator as a function of /3 for projected wave functions with / = 2,6 and a fixed value of 
7 = 50°. The integration has been performed with two sets of integration points in (a, b, c), 
Si = (6, 16, 12) and 5*2 = (16, 16,32). Here, we can observe that for the set S2 the correct 
result of the eigenvalue is obtained for all /3 and J, K. However, the set Si fails both for 
large values of (3 for all /, K and also for smaller deformations with high K = A,Q. The 
poor performance of this choice is clearly seen in the latter case where substantial deviations 
from the correct number are observed. Therefore, as a rule of thumb, the larger the values 
of (/, K) and /3 the more integration points are needed to have good results. The final 
choice will be the one that is able to provide converged results for all (/, K, /3, 7) values. 
Taking into account that the symmetry axis corresponds to pure K = states, one may 
assume that close to the symmetry axis only small i^-components are present. We therefore 
examine the role of the orientation of the intrinsic axes in the PNAMP method. First, we 
explore the convergence of the angular momentum projection using the property given in 
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FIG. 4: (Color online) Expectation values of the total angular momentum operator calculated with 
angular momentum projected states \IK) as a function of the (5 deformation (7 = 50°) and for 
different sets of integration points in the Euler angles {a,b,c) (red circles 52 = (16,16,32), black 
filled boxes Si = (6, 16, 12). The top and bottom panels correspond to / = 2 and / = 6, and their 
corresponding K values, respectively. 

Fig. 3 and projecting symmetric states with the same value of /3 but with 7' = 120° + 7. If 
our assumption is right, we could reduce the number of integration points using instead a 
given wave function an equivalent intrinsic wave function with an orientation closer to the 
K = case. In Fig. 5 we plot as a function of /3 the expectation values of the angular 
momentum operator for intrinsic states with 7 = 50° and also with 7' = 170°. The sets of 
integration points are the same as in Fig. 4. For the set 5*1 with 7 = 50° we observe again 
the loss of convergence whenever /3 and I increase. However, very much improved results are 
obtained for the same set of integration points, 5*1, but projecting the wave functions with 
the 7 = 170° orientation. In addition, the calculation with the set 5*2 reveals the numerical 
origin of the lack of convergence for the set 5*1 with 7 = 50°. Therefore, we will use this 
property to define the mesh in the (/3,7) plane for performing GCM-PNAMP calculations 
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FIG. 5: (Color online) Expectation values of the total angular momentum operator between angular 
momentum projected states / = 0, 2, 4, 6, 8; = and I = 3; K = 2 as a function of the /3 
deformation and for different sets of integration points in the Euler angles and orientation of the 
intrinsic axes, red circles 5*2 = (16, 16, 32), black filled boxes Si = (6, 16, 12) with 7 = 50° and blue 
filled diamonds with 7 = 170°. 

as we will see below. 

The analysis shown in Figs. 4 and 5 has been performed with diagonal matrix elements. 
Since in the GCM calculations we have to consider also non-diagonal matrix elements, we 
have extended our study to this case. We find that in order to ensure a good convergence 
in all cases, the final set of integration points in the Euler angles has to be chosen as 
{Na = 8,Nj, = 16, Nc = 16). We can also exploit the degeneracy illustrated in Fig. 3 
to perform a consistency test of the implementation of the PNAMP method [47]. Using 
symmetry properties of the point group D2 it can be shown, in the notation of eq. 2, that 

{IMK; NZ; = QO°\H\IMK; NZ; = 60°) _ (JOO; NZ; = 180°|i^|/00; NZ; = 180°) 

{IMK; NZ; = 60°\IMK; NZ; = 60°) ~ {m;NZ;^-f = 180° \I 00; N Z; f3-f = 180°) 

(25) 

i.e., the projected energy calculated with a HFB wave function with 7 = 60° is K- 
independent and equal to the projected energy calculated with the HFB wave function 
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K=0 K=2 K=4 K=6 K=8 mixed K=0 

FIG. 6: (Left panel) excitation energies and B(E2) (in e^fm^ units) values for the states before 
and after K mixing -\IMK) and \IM) respectively- with f3 = 0.625 and 7 = 60°. (Right Panel) 
excitation energies and BE(2) values for the state \IMK = 0) with /3 = 0.625 and 7 = 180° 

with 7 = 180°. A similar relation applies for the transition probabilities. In Fig. 6 we show 
the excitation energies and reduced transition probabilities B(E2) calculated with the same 
oblate axially symmetric wave function {/3 = 0.625) but oriented differently in space with 
7 = 60° (left panel) and 7 = 180° (right panel). As expected, we find that the 7 = 60° 
excitation spectrum and transition probabilities are i^-independent and therefore identical 
to the mixed ones. A look to the right panel corroborates also that these quantities co- 
incide with the ones generated with the 7 = 180° intrinsic wave function. Once we have 
analyzed the convergence and consistency of the PNAMP method for a given point in the 
(/3,7) plane we can study the potential energy surfaces (PES) for the different approaches 
(VAP-PN, PNAMP with and without ii'-mixing). We explore first the role of the mesh of 
points needed to cover all different triaxial shapes. Given the better convergence properties 
for wave functions with a large K = component (compare Fig. 5), we divide the calculation 
into two regions, 7 G [0°, 30°] and 7 G [150°, 180°] (see Fig. 7). The last interval is equivalent 
to 7 G [30°, 60°] and we will transform the results to it whenever we plot the different PES 
throughout this paper. Furthermore, the resolution of the PES is affected by the way we 
perform the discretization of the plane. In the lower panels of Fig. 7 we show the VAP-PN 
energy surfaces for a constant step division both in f3 and 7 directions (left part) and for 
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FIG. 7: (Color online) Mesh of points using constant step in /3 and 7 (top left) or triangle division 
(top right) and the corresponding calculated VAP-PN potential energy surfaces (lower panels) 
transformed to the interval 7 e[0°,60°]. The energy is normalized to the minimum of the PES 
(—196.01 MeV) and the contour lines are divided in 1 MeV (black dashed lines) and 2 MeV steps 
(continuous magenta lines). 

a division based on equilateral triangles (right part). The number of points is Npoints = 99 
in both cases. We observe that the distribution of the points in constant steps is not the 
best choice neither for small /3, where for many points almost degenerated states are ob- 
tained, nor for large /3, where a loss of resolution in 7 is observed for increasing values of p. 
It is precisely in this region where the interpolation between distant points produces arti- 
facts or wrong results in the PES such as spurious oscillations, as for example in the region 
(/3 G [1.0, 1.2], 7 e [20°, 40°]) or softening of the contour plots {/3 e [0.6, 1.1], 7 G [50°, 60°]). 
This is rectified with a discretization based on triangles and the results presented hereafter 
are calculated with this mesh. Nevertheless, although only small differences around the min- 
imum of the PES are obtained in the case of ^^Mg, these effects will be enhanced for rather 
7-soft and moderate /3 deformed nuclei. In those cases, the division based on triangles will 
give much better results for the same number of total points included in the calculation and 
will save computing time with respect to the other mesh. 
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V. TRIAXIAL CALCULATIONS FOR ^^MG 



In the previous sections we have studied several aspects needed to ensure a good perfor- 
mance of the full generator coordinate method with the particle number and triaxial angular 
momentum projected wave functions. This previous research is important because the full 
GCM-PNAMP calculation is very demanding in CPU-time and both convergence tests and 
the choice of the relevant parameters should be performed in advance, but nonetheless also 
checked afterwards. In this section the final results for ^^Mg are presented, their calculation 
as mentioned above, have been done with the set of integration points in the Euler angles 
{Na = 8,Nb = 16, Nc = 16). We choose the triangular mesh with Npomts = 99 shown in 
Fig. 7 to solve the constrained particle number projection before the variation (VAP-PN) 
equations. The intrinsic many body wave functions |$(/3, 7)) are expanded in a cartesian har- 
monic oscillator basis and the number of spherical shells included in this basis is Nsheiis = 7 
with an oscillator length of 6 = I.OIA^''^. In Fig. 7 the VAP-PN energy landscape is plotted 
showing a single and well defined minimum at /3 = 0.5,7 = 0° separated by ~ 7.7 MeV 
from the spherical point and ~ 6.1 MeV from the oblate saddle point at /3 = 0.25. These 
results are consistent with the ones obtained in the axial calculation (see Fig. 1) with the 
difference of having a saddle point in the (/3,7) plane instead of a minimum on the oblate 
side. Similar PES are obtained for Skyrme (HFB with particle number projection after 
variation (PN-PAV) included) [47] and Relativistic (BCS without PNP) [48] interactions 
although a softer surface between the spherical point and the minimum is obtained for the 
Skyrme interaction. 



A. Triaxial PNAMP potential energy surfaces and the RVAMPIR approach for 

24Mg 

The solution of the triaxial HWG equation, Eq. 13, does not require to perform a separate 
angular momentum projection in the laboratory system for each component of the GCM 
basis states in the sense of Eq. 22. However, as in the axial case, we expect the PNAMP 
potential energy surfaces to provide insight and a better interpretation of the configuration 
mixing calculations. We can also separate the energy gain due to the triaxial AMP from the 
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FIG. 8: (Color online) PNAMP potential energy surfaces including ii'-mixing in the (/3,7) plane 
for 1 = — 8 and the first eigenvalues in X-space. The PES are normalized to the minimum 
of the surfaces (-200.74, -199.43, -194.04, -196.61, -190.86, -192.27, -186.09, -185.33 MeV for / = 
0, 2, 3, 4, 5, 6, 7, 8, respectively). The contour lines are divided in 1 MeV (black dashed lines) and 2 
MeV steps (continuous magenta lines). States with projected norm less than 10~^ are removed 



one due to the (/3,7) configuration mixing. Furthermore, they are very important because 
the minima of these PES determine the associated RVAMPIR solution. The PNAMP is an 
involved approach that requieres the solution of the HWG equation, Eq. 23, to include the K 
mixing. The HWG eigenstates, Eq. 22, provide real eigenstates of the symmetry operators 
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that can be used, as we shall see below, to generate energy spectra and to calculate transition 
probabilities. 

In Fig. 8 we plot the normalized PNAMP energy landscapes in the (/3,7) plane for the 
lowest eigenvalue in the i^-space for each angular momentum / = Oi" — 81" (see Eq.23). 
In addition, all the points close to the spherical one, and those close to axiality for odd 
values of J, have been removed for J 7^ because their norm is very small. The first 
noticeable aspect is that the VAP-PN axial minimum of Fig. 7 becomes a saddle point, the 
minimum being displaced towards larger /3 values and 7 > 0° for all values of the angular 
momentum, although the barriers between the new minima and the axial prolate saddle 
points are less than 1 MeV. For / = 0i,2i the minima are located in (/3 ~ 0.7,7 ~ 10°) 
while with increasing value of the angular momentum we observe a softening of the PES and 
a displacement of the minimum to larger 7 and smaller /3 deformation, (/3 ~ 0.65,7 ~ 15°) 
for I = Af, 5t and {/3 ~ 0.55, 7 ~ 17°) for / = 6^^, 7^^, 8|. We also note that in the case of 
odd- J values the softening of the PES is in the 7 direction towards the oblate saddle point. 
The energy difference between the VAP-PN and the / = O^" minima is ~ 4.6 MeV while 
the gain in energy due to the inclusion of the triaxial degree of freedom, i.e, the difference 
between the triaxial minimum and the axial saddle point, is ~ 0.7 MeV. Similar results 
have been reported with Skyrme and Relativistic interactions although these studies of the 
PNAMP-PES only extend to / = 0,2 and the effect of increasing triaxiality with growing 
angular momentum has not been analyzed. 

For an interpretation of the configuration mixing calculations it has become customary to 
plot the diagonal matrix elements of the normalized Hamiltonian overlap, see Eq. 14, i.e., 
the JK-projected energy 

^r"(/37) = (26) 

in the (/3,7) plane for the different i^-values. Since the states \IMK; NZ; /S'j) are not 
eigenstates of the angular momentum in the laboratory frame their energies do not have a 
physical meaning. Furthermore the /i^'-projected energy PES and wave functions depend 
on the orientation of the axis in Fig. 3. To illustrate this point we present in Fig. 9 PES's 
calculated in three approaches for different orientations of the nucleus according to Fig. 3. 
We observe that, as expected, one sixth of the circle (for instance 7 = 0° — 60°) is enough 
to describe the PES corresponding to the VAP-PN and the PNAMP ones (corresponding to 
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FIG. 9: (Color online) Potential energy contour plots for ^^Mg in the (/3,7) plane for 7 = 0° — 
360° in different approaches and angular momenta normalized to the corresponding minima. The 
contour lines are divided in 1 MeV (black dashed lines) and 2 MeV steps (continuous magenta lines). 
(Left panels) Top: Particle number projection (VAP), Emin = -196.02 MeV; bottom: PNAMP 
approach for / = Ofi, Emin = —200.74 MeV. (Middle panels) IX-projected energies according to 
Eq. 26. Top: I = 2,K = 0, E^in = -199.42 MeV; bottom: I = 2,K = 2, Emin = -198.78 MeV. 
(Right panels) Lowest eigenvalues of the PNAMP approach. Top: / = 2i, E^in = -199.43 MeV; 
bottom: / = 22, E^in = -195.18 MeV. 

/ = Oi and / = 2i, 22). However for the i^T-projected PES's {I = 2, K = and I = 2, K = 2) 
a semicircle (for instance 7 = 0° — 180°) is needed. Since in the laboratory system all the six 
sectors are equivalent we explicitly see that it is the same to use the region of 7 = 150° — 180° 
than 7 = 30° — 60°, as we have done in the GCM calculations. The contour plots in the IK 
projection can be easily understood looking at Fig. 3. For I = 2, K = 0, the collective AM 
is perpendicular to the z-axis and since semi-classically a rotor will prefer to rotate around 
the axis with the largest moment of inertia it is obvious that the energy minima are around 
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TABLE I: /3 and 7 coordinates of the triaxial PNAMP minima after ii'-mixing as well as excitation 
energies and distribution of K components (i.e., \Hj^^^''^ {f3, 7)^ see Eq. 22 and below) as a function 
of I^. The values of (/3, 7)min may not coincide exactly with those of Fig. 8 because of the finite size 
of the grid used in the calculations. The quoted values are the actual ones used in the /C- mixing 
calculation. The K = ±6, ±8 components, not shown, are exactly zero. 

7 = 0°. For I = 2, K = 2, the collective AM is parallel to the z-axis and in this case the 
minima will be around 7 = 120° and 7 = 240°. Specially for the latter case we see that 
it can be dangerous to make interpretations based on the 7 = 0° — 60° sector. For nuclei 
with more mixing one should also care about the interpretation of the I = 2, K = surface. 
In any case, it is important to note that the K value is not a good quantum number in 
the laboratory frame and therefore it is not an observable. In addition, the distribution 
of K and the corresponding PES can change depending on the orientation of the intrinsic 
wave function (see Fig. 3). Nevertheless, in cases where the i^-mixing is not very large this 
quantum number can be useful to give an interpretation of the different bands that could 
appear in the spectrum. As we will see below, ^^Mg is a very good example of rather pure 
\K\ bands. One should be aware, however, that even with rather pure \K\ = 2 bands, a 
mixing of i^' = 2 and K = —2 takes place and since these states are not orthonormal pitfalls 
may appear. 
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FIG. 10: RVAMPIR spectrum and transition probabilities. 

As discussed in Sect. II B, the minima of the PNAMP potential energy surfaces provide an 
approximation to an angular momentum projection in a variation after projection approach, 
which we have called RVAMPIR. In Table I we present the (/3,7) of the minima of the two 
lowest eigenstates together with the i^'-distribution of the corresponding wave functions. 
As we observe there is almost no mixing, the Oj'", 2^^, 45*", 6^, Sj/" states are K = and the 
2^, SJi*", 42^, 55*", 75*", K = 2. Only at the highest angular momentum we observe very small 
i^'-mixing. This is not the general rule. The amount of K-mixing depends strongly on the 
nucleus and on the 7) point. As mentioned ^^Mg seems to be a nucleus with rather small 
i^'-mixing. The solution of the HWG equation, Eq. 23, at the point (/34in)74in) provides 
the RVAMPIR energies and wave functions of the corresponding states. In Fig. 10 we present 
the energy spectrum and the calculated transition probabilities. Though we will discuss this 
figure in relation with the full GCM results we can compare with the Axial PNAMP-GCM 
excitation spectrum of Fig. 2. The clear difference is the presence of a well developed gamma 
band in the RVAMPIR calculations. 

B. The configuration mixing calculations for ^^Mg 

The final step in the calculation of the spectrum is the GCM-PNAMP method, in which 
simultaneous mixing of the different deformations (/3,7) and K components is performed 
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FIG. 11: (Color online) GCM-PNAMP energies (/ = 2) as a function of the corresponding norm 
eigenvalue, normalized to the highest value, used as cutoff in the definition of the natural basis 
(Eq. 17). 

(see Eq. 1). As we mentioned in Sec. II A, we have to solve the HWG equations separately 
for each value of the angular momentum. These generalized eigenvalue problems are solved 
removing the linear dependence of the states with the definition of the orthonormal natural 
basis (Eq. 17). In order to avoid spurious states in this basis, we have to define a cutoff 
parameter, (, to determine the states in the natural basis, see Eq. 17 and text below. The 
convergence of the triaxial PNAMP-GCM method is showed in Fig. 11 where the lowest 
three energy values obtained for / = 2 are represented as a function of the parameter (. 
Here we distinguish a region of large ( in which the energies are decreasing followed by a 
range of values where the energies are nearly constant. The appearance of these plateaus is 
the signature of the convergence of the GCM method [67]. We observe that this plateau is 
better defined for the 2i and 2^ states as compared to the 2^ . Finally, for small values of C 
the linear dependence shows up and we obtain meaningless values for the energy. The final 
choice for is a value around which we observe a wide plateau for all the levels of interest. 
This value must be kept constant for a given angular momentum in order to guarantee the 
orthogonality of the corresponding wave functions. This analysis has been performed for 
different values of the angular momentum showing in all cases a behavior similar to the one 
presented in Fig. 11. Eventually, we have chosen ( = 10~^ as the final value, similar to the 
one found in Ref. [49]. This procedure can be complemented by an inspection of the shape 
of the collective wave function function of (. 
In the central panel of Fig. 12 we now plot the spectrum of ^^Mg extracted from the 
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FIG. 12: Calculated excitation energies and reduced transition probabilities B(E2) (in e^fm^) 
in ^'^Mg obtained using axially symmetric (left) and triaxial (middle) GCM-PNAMP approaches 
compared to the experimental values (right). The widths of the arrows are proportional to the 
corresponding B values . The experimental values are taken from [68] 

triaxial GCM calculations. We classify the different levels in three bands according to the 
corresponding B(E2) values. The ground state band is formed by a sequence of even angular 
momentum states with a level spacing very similar to that of a rotational band. The first 
excited band consists of states with / = 2, 3, 4, 5, ... as expected for a 7 band. The third band 
is built of even-/ states on top of the second 0^ state. We can also compare the absolute 
value of the ground state energy calculated with different approaches. Evidently, the lowest 
value is obtained with the triaxial GCM-PNAMP method (-201.36 MeV) while -200.74 MeV 
and -200.67 MeV are the results for RVAMPIR and axial GCM-PNAMP approximations 
respectively. Comparing the first two values we observe that the energy gained by mixing 
different shapes is ~0.5 MeV, much less than the correlation provided by PN and/or AM 
restoration. On the other hand, the inclusion of the triaxial degree of freedom within the 
GCM framework gives a similar energy gain (^0.5 MeV) due to the fact that the ground 
state -and also the whole band built on top of it- is already well described by an axial 
calculation in this particular nucleus. Major changes, as we will see below, are however 
found for the excited bands. Concerning the transition probabilities, we observe strong elec- 
tric quadrupole intraband transitions while the interband E2 transitions are much weaker. 
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1.000 _ _ _ _ 
0.922 0.039 _ _ _ 
0.904 0.016 0.032 — — 
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0.834 0.040 0.017 0.016 0.010 


2t 
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ot 
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11.265 
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1.000 — — — — 
0.958 0.021 _ _ _ 
0.795 0.048 0.055 — — 



TABLE II: Decomposition of the norm of GCM-PNAMP collective wave functions into K compo- 
nents (Z)/3,7 \I^K^^''^ if^il)]"^) for the first, second and third bands. The highest values are printed 
in boldface. The excitation energies of the corresponding states are also provided. 

This fact indicates different underlying structures of the bands and the absence of mixing 
between them. We can study the nature of these bands decomposing the collective wave 
functions \F^'^^''^{P,'y)\'^ (Eq. 19) into their K components, 7)^, summing the 

contribution of all deformations (/3,7) for each K. The result is shown in Table II , where 
we clearly observe that the first and third bands are rather pure K = while the second 
band corresponds mainly to li^l = 2 states. Furthermore, we see that for each level the ±K 
components have the same values, as a direct consequence of the time-reversal conservation 
of the intrinsic wave functions. 

The distribution of the states within these bands is supported by the values of the spectro- 
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FIG. 13: (Color online) Spectroscopic quadrupole moments as a function of the angular momentum 
/ calculated for the states in the first (black bullets), second (blue diamonds) and third (magenta 
triangles) bands and for the vibrational-rotational collective model with K = (red boxes) and 
K = 2 (black circles). 



scopic quadrupole moments: 
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where M|^^^ = er^Fj 



are the electrical quadrupole moment operators. In the collective 



rotational model the spectroscopic quadrupole moments for a given |iir|-band take the simple 
form [2]: 



Q,ou{I,K)=Q, 



3K^ - I{I + 1) 



(28) 



(/ + l)(2/ + 3) 

with Qq a constant deformation of the intrinsic macroscopic state. In Fig. 13 we compare 
the triaxial results with the values given for the collective rotational model with K = 
and K = 2 -normalized to J = 2. Here we can clearly observe that the ground state band 
corresponds to a rotational band {K = 0) while the second band matches to a 7-band 
{K = 2) and the third band cannot be described in this simple picture. 
We also plot the probability distribution |F'^'^'^'°"(/3, 7)^ of each GCM state in the (/3,7) 
surface (Fig. 14) summing all the possible K components. The most noticeable aspect is 
that all the states belonging to the same band have a very similar probability distribution in 
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the (/5,7) plane and that the overlap between states of different bands is small. One could 
assume that these facts will lead to the intraband and interband B(E2) values shown in Fig. 
12. However, as we shall see below, the reason for the small interband transitions seems to 
be more related to a i^-hindrance aspect based on the fact that the ground band is a pure 
K = and the 7-band a rather pure ji^'l = 2 band. In particular, all the states in the first 
band have a well defined maximum at {/3 ~ 0.58,7 = 0°) and the probability drops rather 
symmetrically in the /3 and 7 directions. For the second band, the probability distribution 
is concentrated in a region of the plane with {/3 G [0.4 — 1.0], 7 G [0°,35°]) with maxima 
around {(3 ~ 0.7,7 ~ 18°). Finally, the states belonging to the third band show a high 
probability of having spherical shape (02^) or slightly prolate (23", 43", 63") combined with a 
non-negligible mixing of strongly deformed states in the range of /3 G [0.8, 1.3], 7 G [0°, 30°]. 
The PNAMP-PES of Fig. 8 can help to understand the probability distribution of the HWG 
equation. Looking at the ground band panels (O^*", 2^, 4^/", 6^, 8]*") of Fig. 8 we find that all 
show soft triaxial minima close to the axial axis, the contour lines being elongated along the 
radial direction and rather steep along the 7 angle. These states will mix with the mirrored 
ones at 7 = 0°, see Fig. 9, and as a result distributions with a peak at 7 = 0° similar to 
the ones in the left panels of Fig. 14 are expected. If we now concentrate on the panels 
(35*", 5^, 7^), representative of the 7 band, we observe contour lines centered around a soft 
slightly triaxial minimum. These contours, at variance with the ones of the ground band, 
are softer in the 7 angle. We found that the members of the 7 band are rather pure K = 2 
states, that means, the norms of the states along the symmetry axis (7 = 0) are zero. This 
axis acts as a barrier between the states above 7 = 0° and the mirrored ones hindering the 
mechanism described for the ground band. As a result distributions similar to the ones in 
the middle panels of Fig. 14 will be obtained. 

In Fig. 12 we have also compared the triaxial results with axial calculations. In order 
to better understand the results of this comparison, we investigate first the relationship 
between the axial and triaxial collective wave functions. The axial states emerge from the 
7 = 0°— 180° path of the K = component of the corresponding triaxial states. In particular, 
we can relate the ground state bands in both approaches and also the axial 0^, 2^,4^ with 
the triaxial 0^, 2j^, 4^ states (see Figs. 1 and 14). Hence, the comparison between the triaxial 
and axial calculations reveals that both the energies and reduced transition probabilities of 
the ground state band are very similar in both expected. Nevertheless, the small K- 
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FIG. 14: (Color online) GCM-PNAMP collective wave functions {F^'^^^'" {l3,-f)\'^ for the ground 
state (left), second (middle) and third (right) bands, respectively. Contour lines are separated in 
0.01 units. 

mixing for / 7^ lowers the excitation energies of higher angular momentum and therefore, 
the triaxial ground state band is slightly compressed with respect to the axial band. This 
effect, although small, helps to improve the description of the moments of inertia within the 
GCM-PNAMP framework. Larger differences between the axial and triaxial calculations 
appear for the second and third bands. Obviously, the axial calculations are unable to 
describe the 7-band but also the energies and B(E2) for the third triaxial band {K = 
mainly) are modified with respect to the corresponding ones in the axial case. This difference 
is due to both the small i^-mixing and the triaxial configuration around /3 ~ 1.0 that appears 
already for / = (see Fig. 14). 

In Table III we present the average intrinsic deformation parameters and the spectroscopic 
quadrupole moments obtained for the axial and triaxial calculations. In general the average 
P deformation is larger in the triaxial than in the axial calculations. The largest differences 
correspond, obviously, to the states that compose the 7 band in the triaxial case and to the 
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0^ state due to the fake minimum on the oblate side of Fig. 1. Also interesting to notice 
is that though the most probable 7-value in the first band is zero, the average 7 values are 
around 15°. Also the average 7 values for the 7-band are larger than the most probable 
ones. For completeness we also include the values of the spectroscopic quadrupole moments. 
At this point we would like to discuss the bands obtained in the RVAMPIR approach and 
plotted in Fig. 10. At first glance both bands look similar to the corresponding ones of the 
full GCM calculations. A more careful analysis shows that the GCM bands are slightly more 
compressed than the RVAMPIR ones in a better agreement with the experimental results. 
The transition probabilities are also very similar in both approaches. It is really surprising 
that the RVAMPIR is able to provide spectra and transition probabilities comparable to 
the full GCM approach. There are several reasons which explain this behavior. A look to 
Fig. 14 shows that all states of the first and second band, respectively, do have a similar 
probability distribution consisting of one maximum -practically at the same (/3,7) point 
for all /-values- and a homogeneous spread around this point. Such distributions can be 
very well approximated by a delta function at the given point. Furthermore since there is 
no i^-mixing neither in the GCM nor in the RVAMPIR there is no chance that the two 
approaches can differ in this respect. The maxima of band 2 are located practically at the 
same 7) point in both approaches while for band 1 the GCM maximum appears closer to 
the symmetry axis. However since the surfaces are rather flat around these points this does 
not matter too much. With respect to the B{E2) transition probabilities we observe that 
they are rather similar to the ones of the GCM calculations, i.e., they are large for intraband 
and small for interband transitions. Interestingly the 2^, the 2^ and the 3^ RVAMPIR states 
do have the same deformation parameters, see Table I, i.e., we cannot argue, as in the GCM 
case, that the small interband transition probabilities are due to the poor overlap of the 
corresponding wave functions. The reason is that the ground band is a pure K = and 
the 7-band a pure \K\ = 2 band. As a matter of fact, in this case, if we look at Eq. 21, 
we observe that if the factors sandwiched between the collective wave functions do not mix 
strongly the K quantum number, then the transition probabilities are very small. 

Finally we compare the triaxial results with the available experimental data for ^^Mg 
(see Fig. 12). We find a remarkable qualitative agreement between theory and experiment 
both in energies and reduced transition probabilities. In both cases we observe a rotational 
ground state band, a second band associated to a 7-band and a third band with A J = 2. In 
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fact, the theoretical description of the experimentally observed 7-band is one of the major 
achievements of the present model compared to previous implementations. Furthermore, 
in the particular case of ^^Mg, the excitation energies within the ground band are very 
well described even quantitatively with the present calculations as we see in Fig. 12. In 
addition, it is important to emphasize the quality of the theoretical results for the intraband 
and interband reduced transition probabilities which reflect the small K mixing between the 
corresponding bands. Although the improvement of the results with respect to the axial case 
is evident, the band heads of the 7- and, especially, the third band are still calculated too high 
in excitation energy. This is probably due to the lack of the correlations associated to the 
angular momentum restoration before the variation and time-reversal symmetry breaking 
that are not included in this calculation. Additionally, the inclusion of two-quasiparticle 
states would further lower the excitation energies for these band heads. These effects could 
also be present in the ground state bands. However, all these potential improvements are 
beyond the scope of the present work. They would probably lead to a better quantitative 
description of the experimental results although we do not expect qualitative changes in the 
general picture. Research in this direction is in progress. 

VI. SUMMARY 

In summary, we have presented the first implementation of the GCM-PNAMP method 
with fully triaxial intrinsic wave functions found by solving the VAP-PN equations with the 
Gogny interaction. Furthermore, due to the huge computational effort demanded by this 
type of calculations, we have established a protocol for a good performance of the method, 
namely: 

1. Perform first a GCM-PNAMP with only axial K = wave functions in order to choose 
the number of oscillator shells, the relevant interval of /3 deformation and the density 
of points in the collective variable. 

2. Study the convergence of the triaxial angular momentum projection with the number 
of integration points in the Euler angles by looking at the expectation value of P in 
the (/3,7) plane and exploit the symmetries of the intrinsic states. 
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r)spec 
^ax 


1 


0.644 14.81 


0.00 


0.510 





0.00 


2 


0.515 25.98 


0.00 


0.221 





0.00 


2 1 


0.658 13.58 


-20.80 


0.606 





-19.94 


2 2 


0.709 23.16 


21.59 


0.346 


60 


12.70 


2 3 


0.642 16.02 


-19.66 


0.451 





-19.33 


3 1 


0.717 22.05 


-0.18 




4 1 


0.661 13.65 


-26.33 


0.623 





-25.04 


4 2 


0.714 22.60 


-10.96 


0.495 





-20.33 


4 3 


0.541 19.95 


-18.51 








5 1 


0.720 22.67 


-17.51 




6 1 


0.657 14.51 


-28.42 


0.622 





-26.52 


6 2 


0.702 22.92 


-20.97 






-18.08 


6 3 


0.535 17.24 


-23.16 








7 1 


0.730 24.07 


-24.49 




8 1 


0.701 15.36 


-25.96 


0.648 





-28.18 



TABLE III: Average intrinsic deformation parameters and spectroscopic quadrupole moments (in 
e fm) in the triaxial and axial approximations. 

3. Choose a triangular mesh in the (/3,7) plane in order to avoid both redundancy near 
the spherical shape and spurious effects due to a loss of resolution for increasing (3. 

4. Select the converged states as the ones whose energy belongs to a plateau and ensure 
the orthogonality with the other states with the same angular momentum. 

5. Check the convergence of the results in the full triaxial calculation. 

The method has been applied to the study of ^^Mg which has been chosen as a test case 
in previous studies with different interactions. The comparison between axial and triaxial 
results shows minor changes in the ground state band which is predicted to be an axial 
rotational band with X = 0. Only for angular momentum / > 4 some i^-mixing is observed 
giving rise to a small level compression. This result supports the use of axial calculations 
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in these cases. However, the triaxial calculation is also able to reproduce the second band 
associated to a 7-band {K — 2) observed experimentally. 

We have also introduced the RVAMPIR method which provides a more affordable al- 
ternative to the full GCM procedure for the calculation of ground and 7 bands. Wc find 
that this approach provides a good description of the energy levels and the intraband and 
interband transition probabilities for the nucleus ^^Mg. 

Furthermore, the agreement between the theoretical and experimental results is in general 
good although some improvements beyond the scope of this work must be performed in 
order to give a better quantitative description. Some work is in progress in order to take 
into account these effects. In any case, ^^Mg is not a very good example for studying strong 
triaxial effects like the ones mentioned in the introduction and the method will be applied 
in the near future to other systems where both triaxiality and X-mixing play a crucial role 
in describing the experimental data. 
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